library(FactoMineR)
library(factoextra)
library("corrplot")
library(ggsci)
library(plotly)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ComplexUpset)
library(pander)
library(knitr)
library(discoveR)
library(ade4)
library(adegraphics)
library(vegan)
library(MASS)
library(ComplexUpset)
library(RColorBrewer)

Dans tout le jeu de données, on retrouve :
- 9 stress biotiques (Biotrophic bacteria, Fungi, Nématodes, Oomycète, Rhodococcus, Stifenia, Necrotrophic bacteria, Virus et Other biotic).
- 9 stress abiotiques (Heavy metal, UV, Drought, Gamma, Nitrogen, Oxydative stress, Salt, Temperature et Other abiotic).

Seuls les stress biotiques sont étudiés ici. Les fichiers utilisés ont deux colonnes informatives avec le stress appliqué sur l’échantillon en première colonne et le SWAP_ID pour pouvoir remonter aux informations de l’expérience précise.

Le fichier de chaque Gene Set est chargé à partir de la liste fournie dans le script. Les échantillons pour ce set sont séparés enuite pour les stress biotiques et abiotiques.

# Chemin vers dossier fichiers .txt et chargement
setwd("../docs_txt/")
data=read.table('Gene_Swap_NO_NA.dat',header=TRUE,sep='\t')

# Noms des colonnes sans les variables qualitatives
Gene=noquote(colnames(data[,c(-1,-2,-3)]))

#Noms des fichiers pour lecture dans boucle
GOSLIM=c("circadian_rythm.txt","abiotic.txt","biotic.txt","endogenous.txt",
         "external.txt","flower.txt","growth.txt","light.txt","photo.txt",
         "stress.txt")

matrix=cbind(rep(0,17341),Gene) #initialisation matrice

#Remplacer 0 par TRUE ou FAlSE pour chaque gene set
for(i in 1:10){
  datago=read.table(GOSLIM[i],header=TRUE,sep='\t')
  colgo=noquote(colnames(datago[,c(-1,-2,-3)]))
  test=match(Gene,colgo,nomatch =0)
  test[test!=0]<-"TRUE"
  test[test==0]<-"FALSE"
  matrix=cbind(matrix,test) #assemblage des vecteurs pour chaque stress
}


matrix=as.data.frame(matrix)
matrix=matrix[,-1] #enlever première ligne de 0 d'initialisation
names(matrix)=c("Gene","Circadian","Abiotic","Biotic","Endogenous_stimulus",
       "External_Stimulus","Flower","Growth","Light","Photosynthesis",
       "Stress")
goslim=colnames(matrix)[2:11]

#Construction upset
(upset(matrix,goslim, name="GO SLIM", min_size=25,min_degree=1)
+ ggtitle("Tous les gènes (173441)"))

setwd("../docs_txt/")

#Variables importantes pour code, noms des fichiers pour lecture, nom des gene set pour entêtes personnalisées, différents stress..
files=c("circadian_rythm.txt","abiotic.txt","biotic.txt","endogenous.txt",
         "external.txt","flower.txt","growth.txt","light.txt","photo.txt",
         "stress.txt","random.txt")

name=c("Circadian rythm","Abiotic","Biotic","Endogenous stimulus",
       "External stimulus","Flower","Growth","Light","Photosynthesis",
       "Stress","Random")

abiotic=c("DROUGHT","GAMMA","HEAVY.METAL","NITROGEN","OTHER-ABIOTIC","OXYDATIVE.STRESS","SALT","TEMPERATURE","UV")

#Boucle principale lecture Gene set
for(i in 1:11){
#Chargement Gene set 
  data=read.table(files[i],header=TRUE,sep='\t')
  #Conservation des échantillons seulement abiotiques
  data_abiotic=data[data[,2] %in% abiotic, ]
  row.names(data_abiotic) <- data_abiotic$Swap_ID
  
  #ACP avec données abiotiques
  res1=PCA(data_abiotic[,c(-1,-3)],scale.unit=FALSE,graph=FALSE,quali.sup=1)
  fic.pca1<-dudi.pca(data_abiotic[,c(-1,-2,-3)],center=FALSE,scale=FALSE,scannf=FALSE)
  
  #BCA données abiotiques
  fic.bca1<-bca(fic.pca1,fac=as.factor(data_abiotic$vec),scannf=FALSE)
  fic.tst<-randtest(fic.bca1)
  
#--------------------------------------------------------------------------------------------------#
                                        #Ecriture dans ouput HTML
#--------------------------------------------------------------------------------------------------#
  
  cat("\n")
  pander::pandoc.header(paste0("**Gene set : ",name[i]," (",dim(data)[2]-3," gènes)**\n"), level = 1)
  
#--------------------------------------------------------------------------------------------------#
                                                   #ACP
#--------------------------------------------------------------------------------------------------#  
  cat("\n")
  pander::pandoc.header(paste0("**ACP (SWAP_ID)**\n"), level = 2)
  
      g1<-fviz_pca_ind(res1,
             repel=TRUE,
             habillage = 1, # color by groups
             addEllipses = TRUE, # Concentration ellipses
             palette="Set1")
      print(g1)
      
      
      cat("\n")#Screeplot intertie
      pander::pandoc.header(paste0("*Inertie*\n"), level = 3)
      print(fviz_screeplot(res1, addlabels = TRUE, ylim = c(0, 50)))
      
      
      cat("\n")#Screeplot contribution individus et variables
      pander::pandoc.header(paste0("*Contribution*\n"), level =3)
      
      cat("\n")
      pander::pandoc.header(paste0("**Individus**\n"), level =4)
      print((fviz_contrib(res1, choice = "ind", axes = 1, top = 10)
        +
      fviz_contrib(res1, choice = "ind", axes = 2, top = 10)))
      
      cat("\n")
      pander::pandoc.header(paste0("**Variables**\n"), level =4)
      print((fviz_contrib(res1, choice = "var", axes = 1, top = 10)
        +
      fviz_contrib(res1, choice = "var", axes = 2, top = 10)))
      
      
      cat("\n")#Boxplot des meilleurs contributeurs 
      pander::pandoc.header(paste0("**Distribution meilleurs contributeurs**\n"), level =4)
      #Récupération des labels des contributeurs via les plots factoextra
      #Découpage en 2 fois pour les deux dimensions
      p=fviz_contrib(res1, choice = "var", axes = 1, top = 10)
      p2=fviz_contrib(res1, choice = "var", axes = 2, top = 10)
      contrib=p$data[order(p$data$contrib,decreasing=TRUE), ]
      contrib2=p2$data[order(p2$data$contrib,decreasing=TRUE), ]
      contrib=contrib[1:5,]
      contrib2=contrib2[1:5,]
      
      #Match entre noms des colonnes data_abiotic et les noms des meilleurs contributeurs
      col_names=colnames(data_abiotic)
      index=match(contrib$name,col_names)
      index=index-2
      index2=match(contrib2$name,col_names)
      index2=index2-2
      
      #Récupération des données de puces (sans colonnes qualitatives) pour les colonnes d'intérêt
      contrib_data=data_abiotic[,c(-1,-2)]
      contrib_data=contrib_data[,index]
      contrib_data2=data_abiotic[,c(-1,-2)]
      contrib_data2=contrib_data2[,index2]

      
      data_long=contrib_data%>%pivot_longer(
                       cols = c(1:5),
                       names_to = "Gene",
                       values_to = "Expression",)
      
      g1<-ggplot(data_long)+
          geom_boxplot(aes(x=Gene, y = Expression))+
          theme(axis.text.x = element_text(angle = 60, hjust = 1))+
          ggtitle("Contributions Dim 1")
  
      
      
      data_long2=contrib_data2%>%pivot_longer(
                       cols = c(1:5),
                       names_to = "Gene",
                       values_to = "Expression")
      
      g2<-ggplot(data_long2)+
          geom_boxplot(aes(x=Gene, y = Expression))+
          theme(axis.text.x = element_text(angle = 60, hjust = 1))+
          ggtitle("Contributions Dim 2")
  
      #Détermination des limites en fonction des valeurs des deux dimensions pour rendre comparaison plus facile
      g1<-g1+ylim(min(data_long$Expression,data_long2$Expression),
                  max(data_long$Expression,data_long2$Expression))
      g2<-g2+ylim(min(data_long$Expression,data_long2$Expression),
                  max(data_long$Expression,data_long2$Expression))
  
      print(g1+g2)

#--------------------------------------------------------------------------------------------------#
                                                #BCA
#--------------------------------------------------------------------------------------------------#
      
      cat("\n")
      pander::pandoc.header(paste0("**BCA**\n"), level = 2)
      cat("\n")
      pander::pandoc.header(paste0("*Toutes positions*\n"), level = 3)
      cat("\n")
      
      .z<-paste("BCA sur position \n ratio= ",round(fic.bca1$ratio,2)," pval= ",signif(fic.tst$pvalue,digits=2))

      s.class(fic.bca1$ls,
              fac=as.factor(data_abiotic$vec),
              col=brewer.pal(n=9,name="Set1"),
              psub=list(text=.z,cex=1,position="topleft"),
              plabels.cex=0,
              key.cex=0.8,
              key.text.cex=0.8)

      cat("\n")
      pander::pandoc.header(paste0("*Par position*\n"), level = 3)
      
      #Plot pour chaque stress
      plotlist=list()
      for(j in 1:9){
        s.class(fic.bca1$ls,
                fac=as.factor(data_abiotic$vec),
                psub=list(text=paste("BCA \n", abiotic[j]),cex=1,position="topleft"),
                plabels.cex=0,
                key.cex=0.4,
                key.text.cex=0.4,
                col=c(rep(grey(0.95),j-1),brewer.pal(n=9,name="Set1")[j],rep(grey(0.95),9-j)),
                plot=FALSE,
                chullSize=0,
                ellipseSize=0)->plotlist[j]
      }
      
      ADEgS(plotlist,layout=c(3,3))
}

Gene set : Circadian rythm (57 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Abiotic (661 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Biotic (420 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Endogenous stimulus (523 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : External stimulus (552 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Flower (186 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Growth (222 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Light (292 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Photosynthesis (75 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Stress (1177 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Random (50 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position